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Figure Captions 


Fig. 1 Location of Magnetic Survey Observations available 1955-1962. 

Only one point plotted in each 0.5° x 0.5° latitude - longitude 
block. 

Fig. 2 Location of data selected 1940-1962 for field analysis. Shaded 
areas designate those for which some data are available for 
the years 1900-1939. 

Fig. 3 (a-d) Contours of the geomagnetic field in gauss (F, H, Z) or 
degrees (I), synthesized from the coefficients in Table 4 for 
epoch 1965 . 0 . 

(a) F (Total Intensity). All centers are 'Highs' except the South 
American 'Low' of 0.238 T. 

(b) H (Horizontal Intensity). The + signs are the dip pole 

positions (north: 75*6° N, 101°W; south: 66.3° S, 141° E) . 

The two centers near the equator are 'High', the center near 
the southern tip of Africa is a 'Low' . 

(c) I (Inclination) 

(d) Z (Vertical Intensity) 

Fig. 4 (a-b) Contours of the secular change of the geomagnetic field in 
gamma/year synthesized from the coefficients given in Table 4 . 

(a) F (Total Intensity) secular change 

(b) Z (Vertical Intensity) secular change 

Fig. 5 Magnetic survey tracks selected for comparison with computed 
fields. Shown are data positions from the U. S. aircraft 
project MAGNET, the U. S. ships Rehoboth and Verna , the Japanese 
ship Soya , and the Russian ship Zarya . 
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Fig. 6 (a-d) 


Fig. 7 (a-c) 


(a) 

(b) 

(c) 


Plots of measured minus computed values of the magnetic 
field components taken over the tracks illustrated in 
Fig. 5« The ordinate scale is marked in 400y increments. 
Model fields used are A from Table 4 and LME as referenced 
in Table 7. (no separate captions) 

Comparison between observatory annual means (dots) and 
field components computed using coefficients A (solid 
lines) and LME (dashed lines) . 

Sodankyla, Finland: 67° 22.0' N, 26 0 39*0' E (1914-1944); 

67° 22.1' N, 26° 37.8' E (1946-1961). 

San Juan, Puerto Rico: 18° 22.9' N, 66° 7.1' W (1926-1963) 

Alibag, India: 18° 38.3' N, 72° 52.3' E (1904-1961). 


Abstract 


A new determination is made of the geomagnetic field and its 
secular change using a set of 21695 selected survey and observatory 
annual mean data available for the interval 1940-1963* The field 
is given by a series of 63 spherical harmonics ( g ® and h™ to n=m=7) 
and 35 time derivatives (g^ and h™ to n=m=5)* Although it is inferred 
from these results that a better fit should be possible with existing 
data if more coefficients are used, the improvement is sufficiently 
substantial to recommend that it be used as a replacement for the 
previously derived coefficients ( Jensen and Cain , 1962) . One re- 
finement in the present determination includes the use of the oblate 
earth in place of the spherical approximation. It is illustrated 
that this modification is necessary to obtain meaningful extrapolations 
above the earth's surface. 

The standard errors and confidence levels are given for the coef- 
ficients lo show that almost all are well determined. The rms devi- 


ations of this fit from several selections of survey data are used in 
comparison with similar results from other coefficient sets to show 
that this set of coefficients is considerably more accurate than any 
published previously and is comparable to a recent determination by 
Leaton et al (1964) • 




Introduction : The history of the determination of the main geomagnetic 

field goes back over a century with excellent bibliographies to be 
found in such works as those by Chapman and Bartels (1940a) > Vestine 
(1947), Mauersberger (1952). Fanselau (1956), and Kautzleben (1963). 

Some of the more recent publications on the subject include those by 
Vestine (I960), Euler (1963), Fougere (1963a), Heppner (1963), Leaton 
(1963) , Vestine et al (1963a, b) . 

In this long history the techniques of evaluation of the main field 
have generally included first determining magnetic charts from the data 
and then performing a spherical harmonic analysis on equi-spaced grid 
points from the charts. The basic evaluation of the data has thus 
occurred in the map-making phase of the processes, rather than in the 
spherical harmonic analysis itself. These maps have generally been 
hand-contoured from a collection of magnetic survey data that have been 
extrapolated to the epoch of the map by attempting to interpolate the 
secular change in the field using the set of data from magnetic obser- 
vatories and repeat stations. Although the data analysis is made very 
difficult by the poor distribution over the earth of the total set of 
survey and secular change data, the accuracy of the final results has 
also been limited by the errors introduced in producing the charts. 

It was with the idea of instead determining a fit directly to the 
survey data themselves that Jensen and Cain (1962) made a computation 
of the main field using a selection from 74000 H and F data available 
for the period 1940-1960. Although the resulting set of spherical har- 
monics was an improvement over the previously available coefficients 
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(cf. Leonard, 1963) for the field near its epoch of I960, there were 
still serious discrepancies between the computed field and measure- 
ments made near thp earth's surface. . 

Using similar techniques of curve fitting, a set of spherical 
harmonics was also computed (Cain, et al, 1962 ) from the Vanguard 3 
magnetic data. As pointed out by Cain and Hendricks (1964) the mere 
fitting of magnetic field data covering a given volume with a set of 
spherical harmonic coefficients does not necessarily imply that the 
main field is correctly evaluated elsewhere. 

Assuming that one wishes to compute a set of spherical harmonic 
coefficients representing the main field directly from the total set 
of available data there are several approaches that can be taken. The 
analysis can either include the determination of the coefficients of 
secular change along with the spatial terms, or instead, the data can 
be adjusted to a particular epoch by a separate analysis of secular 
change and then the spatial analysis made on the adjusted data. The 
only data of sufficient quality to attempt a determination of secular 
change are the annual mean values from the set of about 100 magnetic 
observatories. However, since these observatories are mainly confined 
to the major land masses and are predominently in the northern hemis- 
phere, a determination of secular change from them alone cannot be 
representative of secular change over the whole earth. The quality of 
repeat station data is sufficiently low, particularly in remote areas, 
that it is unlikely that their inclusion would aid substantially in 
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improving this determination. For this reason the approach followed 
here has been to perform the analysis on uncorrected data making the 
simultaneous determination of the time derivatives of the spherical 
harmonic coefficients along with the spatial terms themselves. In 
this way both the observatory and survey data contribute to the de- 
termination of secular change. 

A second basic decision to be made in the determination of the 
main field is whether to use X, Y, Z as derived from the original ob- 
servations of such components as D, H, I and F or to use the observations 
themselves. This is a difficult choice since the use of X, Y, Z data 
requires only the solution to linear equations whereas the use of 
angular data D and I or combined component data H and F requires that 
a non-linear analysis must be performed. However, with the present 
distribution of data where the observations in some of the remote 
regions are only of F it is clear that it is not always feasible to 
compute the vectors X, Y and Z. Also, the only presently available 
satellite data is F (Vanguard 3) which should be used since they are 
the only data which have an appreciable distribution in altitude. A 
further problem with the use of X, Y and Z data is that not only are 
the errors more difficult to estimate, but also they can be systemati- 
cally enhanced over different parts of the globe. For example, Z may 
be obtained from observations of the horizontal intensity H and the 
inclination I, from Z=H tan I. The inequality |6ZU;6H tan lj + }Hsec 2 161 j 
shows that the errors in the computed Z are commensurate with those of 
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I in low latitudes and reflect those of both H and I in middle latitudes. 
However, writing the second term as j(F6l)/cosl| it is apparent that the 
errors in Z reach large proportions in high latitudes as cosl-0. 

Thus the decision was made here to use only observed components in 
the analysis and also to attempt weighting them according to their 
estimated accuracy. This choice implies that if all three components 
are measured, they are treated as three separate observations and in- 
cluded independently into the fit. The weighted sums of squares of the 
quantities (measured minus computed) AX, AY, AZ, AH, and AF are minimized. 
The angular quantities AD and A I are taken into account by converting 
to force 'units using the additional weighting factors H and F (computed) 
respectively. 

A third decision in the analysis of the main field is whether to 
add the refinement of including the earth's oblateness. As we will 
illustrate subsequently there are compelling reasons to include the 
earth's shape if one expects to realistically extrapolate the result 
above the surface. We thus chose to make this deviation from past 
practice even though the resulting coefficients will not be directly 
comparable with those determined previously. There will be a slight 
additional complication in their use since one must take the earth's 
shape into account for synthesis of the field even at the earth's surface. 
2. Method 

We here follow the usual spherical harmonic potential expansion 
(cf. Chapman and Bar tels . 1940 p. 639) 

v = a J 1 { (J) n T n e + (|) n+1 } 



- 5 - 


where • T n = | Q (gjj cos rap + h™ sin rap) p“ (9) , 
r, 9 ,cp = spherical coordinates corresponding to geocentric radius, 
colatitude, and longitude respectively, g™, h® = Gauss (Schmidt 
normalized) coefficients (g° < 0) , and P® (9) = Schmidt's quasi- 
normalized polynomials. The field can be derived by a straight- 
forward evaluation of F = - VV. Although the T® terms are included 
here in the expression for V, they were neglected in the determination 


ws j-ll uc; cl* 


-his omission was deliberate 


for a number of practical and theoretical reasons. First, it is now 
estimated that the true value of such external terms could only con- 
tribute at the level of a few tens of gammas at most from either the 
ring current (Akasofu, Cain and Chapman . 1962 j Hoffman and Bracken . 
1964) or the cavity field ( Mead . 1964, Midgley . 1964) and would thus 
be very difficult to determine in the anomaly field 'noise' in the 
data. Secondly, as will be shown subsequently, the computer core 
size was already exceeded by the number of significant coefficients 
required to fit the internal field. 

Now let one component of the field (X, Y, Z, D, I, H or F) be re- 
presented by the functional form: 

C = C(g, r, 9, 9 , t) 

where g is now the set of all coefficients (e.g. g®, h™, g™...) and r, 
9, 9 and t are the coordinates of the observation point. If there is 
given a set of observations C^ we then attempt to find a set of g's to 
make the sum of squares of the difference between the observed Cj. and 
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computed C (at the cpordinates of the observation) a minimum. The 

quantity co be minimized x 2 is thus: 

X 2 = E CC ± - C(g, r., 0., cp i , t.] 2 
i 

where the summation on i is taken over all component observations. 
Following the usual least squares method we differentiate in turn by 
each of the g's and solve the resulting expressions simultaneously 

ix 2 = 0 = -2? [C i -c] \ k = 1, 2, ...n 

g k 1 ^ ^ g kJ 

However, if the function C is obtained from the spherical harmonic 

expansion of the geomagnetic potential function, the above set of 

equations is non-linear in the g's, and direct methods of solution 

are unknown if not impossible. A method of successive approximation 

was formulated to avoid this problem. 

Assuming that some approximate values of the g's are available, the 

expression for C can be expanded into a Taylor's series to first order: 

~ n 

C = C(g, r,9, cp, t) = C (g, r, 9, cp t) + Z 

k = 1 

where C and ( \ are the function C and its derivative with respect 

l * g Vo 

to g evaluated using the approximate values of the g's, and 6g^ is the 
correction to be applied to the g^ parameter. If the least squares 
procedure is now followed, the resulting set of equations is linear in 
the Sg^'s can be readily solved. 

Before writing the normal equations we first modify the expression 
for x 2 to include the weights for different types of observations and 
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different components as follows: 

X 2 = S [C. - c(g, r., 0., cp., t.)] 3 w. 
where the weights are applied to each component observation as 
w=l/a and a is the estimated accuracy (in gammas) for that class of 
observation. For the components D and I these estimated accuracies 
(in radians) are multiplied by (the computed) H and F respectively. 


Cn substituting for C, 

X s = ? [Ci-C.-E 


Ug, 


0 u 




_r.e .ease scuares croceaure is tc minimize witn resce; 


:t tc 6g k - 


Differentiating with respect to 6g k and setting the partial derivatives 
equal to zero gives 


n 


JUS? = 0 = -2 E [C,-0 -E fj_C \ fig,] /J_C 'l 

M6g k ) 1 J=i 


for all k. 

The resulting kth normal equation is then 

E fig. E w. f^C \ (JlS. \ = E w. (C.-i 

ii 0 i 1 lf«, A TO. i 1 1 


O (te. 
lit 


k / o 


If we now let 



= Sw i 
1 



w k = E "i t°i - C o> 


and 



- 8 - 


we may write the normal equations as 
n 


£ 

j=l 


V 



from which 

6g * = |i w * 

where D” 1 are the elements of the matrix inverse of D. 

Of course such a computation is not meaningful unless some measure 
of the error of the results is evaluated. If we assume that the errors 
in the measurements C i are normally distributed, we can draw some con- 
clusions about the coefficients found. 

We know that a 2 = {£ (AC^ ) 2 Wj_} /£ w- 

i i 


is a minimum variance estimator of the true variance. Secondly, 

^ - s§V g . 

has a "t" distribution where a?. = o , r and g. is one of the co- 

& i 11 ± 

* m 

efficients g£, h®, g®, . Thus if we wish to test the hypothesis 

that g d 4 0, we must set g? = 0 and compare g.j/a d with the value of 
the "t" distribution for any given confidence level. For the present 
case where the number of degrees of freedom, N data points “N coe ffi c ients 
is large the "t" distribution is given to adequate accuracy by its 
asymptotic form, the normal distribution, so that t^ = 1.96 and t^ 0 =0.67 
are the 95$ and 50$ confidence levels respectively. If g.j/a gi ^ *95, 
then with 95$ confidence / 0. Each g^ can thus be said to lie within 
the range + 1.96cr g . with 95$ confidence and the narrower range + 0.67 a gi 
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with 50% confidence. Although these estimates are statistically valid 
they must be viewed with some additional information to decide whether 
a particular spherical harmonic makes a meaningful contribution to the 



separate components of a vector, we have chosen to estimate the sigma 
of the magnitude R™ = 
for g and h. In addition, the set of R/cr^ is considered as a whole for 
a given degree n before assessing its contribution to the expansion. 

A further comparison is made between the o's and the comparative values 
of g k using different selections of data. A truncation of the potential 
expansion is thus made only on n based on these two sources of information. 

The formulation up to the derivation of the is basically the 
same as that used by J ensen and Cain (1962) . The refinements added in the 
present work include the addition of the weights w^ and the extension 
of the procedure to include components other than H and F. Also, in the 
prior analysis there were some initial difficulties with the computation 
and inversion of Dj k that were erroneously thought to be due to the 
size of the matrix. In this previous "calculation instead of determining 

the g k 's for all values of k (i.e. for g°, g^, h^, g^ gn, b^, where 

n = m = maximum degree and order of the spherical harmonic expansion) , 
an iterative procedure was used whereby sets of only 16 or less 6g k 's 
were computed at once holding the rest constant. That is, corrections 

were made first on g^, g^, h^, gy h^ holding all other g's and h's 

constant. Another pass through the data was then made correcting only 


\/g m2 4 


as the rms a of the individual ct 1 s 
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o 1 12 4, 

on the n=4 (g. , g , — h ,h — h ) terms. Each degree harmonic was corrected 

4 4 4 4 4 

in turn up to the maximum and then the process repeated until the 
corrections became small. Although improvements in existing sets of 
coefficients were in fact deduced there was sufficient correlation 


between the various harmonics to cause the resultant solutions to never 


quite stabilize. For the present calculation, however we adopted the 
procedure of solving for all the values of k simultaneously. This tech- 
nique was found to have good results if enough data points were used and 
if the matrix of values entering the normal equations was formed and 
solved in double precision (17 digit) arithmetic. The number of data 
necessary appears to depend jointly on the distribution and type of com- 
ponent and the accuracy of the initial estimates for g^. F data are 
troublesome if the original estimate is bad or if the data are poorly 
distributed. Component data do not need to be well distributed. The 
6g k become small (~10“^g k ) after only two passes through the data 
whereas using the prior procedure the g k continued to change slowly and 
systematically for the correlated coefficients even after several iter- 
ations . 

A further refinement in the present computation scheme over that 
used by Jensen and Cain (1962) is the allowance for using a spheroidal 
earth instead of assuming sphericity. This is done by entering into 
the normal equations for each datum the geocentric coordinates r, 9 
instead of the geodetic colatitude and an r derived from the mean earth 
radius. Most previous works have used 9 = 90°-\ and r = 6371. 2+h where 
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X is the geodetic latitude and h the altitude of the observation above 
the earth's surface. We instead now use the relations 

cot 0 = f h VA^cos 2 X+ if sin 2 X +Ep ] tan X and 

h «/A 2 cos 2 X+]fsinX +A 2 

r^ = h 2 + 2h */a 2 cos 2 X+fiP sin 2 X +A 4 cos 2 X+B t sin 2 X 

A 2 cos 2 X+eP sin 2 X 

for determining 0, r given h and X, where A and B are respectively the 
earth's equatorial and polar radius. The value used in the potential 
expansion V for a was 6371.2 Km to correspond to the earth's mean radius. 
Thus in recomputing the field from the derived set of coefficients a 
truly geocentric r and 9 must be used. The correction for r is the more 
significant of these whereas the correction for the latitude is quite 
small. To be exactly consistent with the directions of the measured 
data on the spheroidal earth (e.g. X and Z) a small rotation also needs 
to be made from the geocentric directions derived from F = -VV. That 
is, if the geocentric components of F are given by 

F r = zAl, F q = -I Al> ^d F = -_1 i_V, 

,X r r 3 9 ^ r sin9 b <p 

then the components X, Y and Z are given by 

X = -F q cos6 -F r sin6 



Z = Fq sin6 -F r cos6 

where 6 is the angle between geodetic and geocentric latitude and can 
be computed from expressions such as sin6 = sinX sin9 - cosX cos9. Since 
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6 is only of the order of 0.2° the correction is not of much consequence 
for most evaluations of the field at the earth's surface. 

Both Schmidt ( Chapman and Bartels . 1940 , pg« 641) Jones and 
Mel o tie (1953) contended that the use of a spheroidal earth in these 
calculations was not a meaningful improvement. However, the accuracy 
of these past workers was only of the order of one percent. Further 
they did not need to extrapolate the field to any significant distance 
above the earth's surface. Making the approximations r = 6371.2 
+ h and 9 = 90° - X during the analysis for g^ and h™ meant that in 
the synthesis of the field components one also used the same relations. 
Tlu.s is mathematically equivalent to mapping the field from the spheroid 
onto the mean sphere in finding the coefficients and then mapping from 
the sphere back to the spheroid for their evaluation. The inherent 
errors in this process can best be illustrated by the following numerical 
example. 

A set of spherical harmonic coefficients was taken to represent the 
exact potential function, V (r, 0, cp). The form of this potential is 
not important other than that it used 48 coefficients and did represent 
the actual geomagnetic field to about 2 %. Using this potential, a set of 
"observed" data was calculated to simulate exactly what an observer 
would measure if he were to occupy each point with his magnetometers, 
levels, and surveying instruments on the spheroidal earth. The "observed" 
F, H, X, Y, and Z were computed over a 10° mesh in longitude and latitude, 
and at altitudes of 0, 10, 1000, 10,000 km. Using the set of Z's for 
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zero altitude, a set of coefficients was calculated assuming a spherical 
earth exactly as has been done by many investigations in the past. In 
order to make this fitted field agree with the "observed" field to less 
than one gamma error, it was necessary to compute g^'s and h® ' s up to 
n = 16, m = 16. The difficult areas to fit were near the poles. 

Now using this set of coefficients, "computed" values of the field 


were calculated for 

the same points 

as the 

“observed" set, now assuming 

a spherical earth, 
the largest errors 

These two sets 
are given in the 
Maximum errors 

of field 
table : 
kfifl ma) 

points 

were compared and 

Altitude 
(km. ) 

Z 

X 

Y 

H 

F 

0 

0 

154 

25 

153 

149 

10 

0 

152 

25 

151 

148 

1000 

49 

87 

13 

86 

65 

10,000 

10 

7 

1 

7 

10 


When the absolute magnitudes of the field are taken into account 
the result is that the maximum error at all altitudes is of the order 
of 0 .%. This percentage error would be the approximate inconsistency 
that would constitute a 'noise' level if the shape of the earth were 
ignored in a field determination. 

3 • Selection of Data 

Another decision necessary in computations such as this is whether 
to perform some reduction and smoothing of the data prior to analysis 
or instead to analyze the raw data. As previously dicussed most of the 
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prior determinations of the main field have been done using compiled 
charts. This technique represents the ultimate degree of data reduction 
which utilizes little of the potential form of the field and raises 
many questions which are more properly left to the data analysis. For 
example, map data have frequently been generated for the earth's sur- 
face by merely using a dipole relation to "reduce" airborne observations. 
Such procedures obviously destroy what little r dependent data are 
available for aiding in the determination of the potential terms. 

There are, however, valid reasons for considering some reduction 
of the survey data prior to analysis. The main theoretical incentive 
is that in making a spherical harmonic analysis determining a non- 
infinite set of coefficients from an irregularly distributed set of 
data there is a tendency for the neglected higher order harmonics to 
influence the values of the lower order terms. This 'aliasing' ( Blackman 
and Tukey , 1958 ) might be reduced by spatially smoothing the data so 
that only those wavelengths are left which correspond to the degrees of 
harmonics n determined. Thus in any analysis of a given set of data 
the values of the lower order terms change somewhat as more coefficients 
are added to the series unless such smoothing is done. However, the 
influence of the high order irregularities can be minimized, particularly 
for the lowest degree terms, by carrying the determination to a suffi- 
ciently high degree. This procedure is practicable so long as the 
neglected harmonics are of sufficiently low amplitude. Fortunately, 
as shown by Alldredge ( 1963 ) the amplitudes do in fact appear to decrease 
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at least for degrees of n continuing beyond 10. 

There are, of course, practical reasons for wishing to reduce the 
set of data lo be analyzed, in that the computation time would be 
impractically large if a numerical fit were made to all available 
survey data in one computation. However, the nature of the data them- 
selves make the decisions as to the type of smoothing very complex. 

For this reason we have here chosen to instead make a selection from 
the currently available data with simple criteria based on the result 
desired. Since the original intent of the study was to develop the 
most accurate reference field for the current epoch for use in con- 
junction with the analysis of data from satellites to be launched during 
the IQSY, the l atest set of survey observations that give an adequate 
coverage of the earth was selected. 

The set of all data available from the Geomagnetism Division of 
the U. S. Coast and Geodetic Survey for the period 1900-1962 comprised 
215,757 observations of one or more components or a total of about 
450,000 component observations ( Hendricks and Cain, 1963). The per- 
centage distribution of these data by decade is given in Table 1 and by 
altitude in Table 2. The main bulk of the latest data are from air- 
borne observations ( Serson , 1957, I960; USN00 . 1963; Behrendt and Wold , 
1963). The only satellite data included in this set are from Vanguard 
3 (Cain, et al 1962) . A plot of the positions of all data taken since 
1955 is given in Figure 1. Only one observation per 0.5° by 0.5° block 
of longitude and latitude is represented on this diagram. As i„- ^vious 


Table 1 


Distribution of magnetic survey 

Altitude (Km) 

0 (surface) 

0.1-2 (airborne) 

2-4 (airborne) 

4-8 (airborne) 

510-3750 (Vanguard 3) 


Distribution of magnetic survey 
Y ears 

1900 - 1909 
1910 - 1919 
1920 - 1929 
1930 - 1939 
1940 - 1949 
1950 - 1959 


data (1900-1962) by altitude. 

Percent 

55 

4 

32 

8 

1 

Table 2 

data by Decade 

Percent 

7 

12 

7 

6 

25 

27 


1960 - 1962 


15 




F igure 
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from this Figure any determination of the earth's field which would 
be representative of that over the whole earth would require data 
to fill in the large gaps over the southern oceans and Asia. 

In order to provide a reasonably uniform distribution of data 
and to reduce the number to that practical for analysis, an equal-area 
grid was constructed based on the area of a 2° x 2° latitude -longitude 
block at the equator. The earth was thus divided into 180 lunes of 
longitude and segments of latitude within each lune so the area would 
be the same as that for the lowest latitude blocks bounded by the 
equator and the + 2° parallels. There were thus 58 segments in lati- 
tude from pole-to-pole with the longest of almost 15° near the poles. 

The latest survey or annual mean observations were then selected from 
all data for each such block. It was surprising to find only 8112 of 
the total 10440 such blocks filled by data after 1900! A time-area 
study of the resulting selected data was then made and it was found 
that a relatively small quantity would be lost by deleting all data 
prior to 1940. The main reason for this fact is that the latest ob- 
served data available for the Asian continent is given for epoch 1940.0. 
The only other substantial area that would be lost by restricting the 
data to post 1940 observations was the South Pacific. However, these 
South Pacific data are now almost worthless since they were taken over 
forty years ago and there are no magnetic observatories in adjacent 
areas to aid in relating them to present data. 
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One danger in making a fit to the resulting data set was that it 
would be difficult to distinguish secular changes from the ordinary 
spatial variations. In order to connect the different areas in time 
a selection of observatory annual means was added from the period 
1940-1962. To avoid distributing these too densely in areas of high 
observatory density the number of annual means was limited to one per 
year per 10° by 10° block of latitude and longitude. 

Figure 2 is a plot of the positions of the data actually used in 
this computation. The shaded areas indicate the largest of those for 
which data were available only prior to 1940 and were hence discarded. 

In preparing these data for analysis an effort was made to select 
from each data source only the observed components. Since each compon- 
ent is fitted separately the use of derived data would have placed too 
high a weight on a given observation. In instances where no choice was 
obvious from the information at hand a maximum of three components was 
selected. 

The set of weights w^ previously defined was then assigned to various 
data classes according to the estimated accuracies of the observations. 
This assignment of weights as well as the selection of the observed 
components was made with the cooperation of the Geomagnetism Division 
of the U. S. Coast and Geodetic Survey. These estimates cannot be 
considered more than qualitative since they are complicated by many 
factors. For example, in the previous formulation there is no allowance 
for errors in the position (or time) for a given observation. Thus the 




Fig. 
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weights must include some allowance for positional errors. This factor 
is most clearly shown by the Vanguard 3 data where it is noted by Cain 
et al (1962) that whereas the errors in the actual observation were 
only of the order of a gamma, the several kilometer error in position 
provided an inherent 'noise' level of the order of lOy to the data. 
Positional errors become even more important near the earth's surface 
where the spatial gradients are larger. The gradients that are signif- 
icant are of course not those of the crustal anomalies but those of the 
derived field. The weight assignments also depend on the smoothing 
involved in the observations. This smoothing is either done spatially 
by averaging observations over a track segment of the observing vehicle 

'*'W >~>w' 

(e.g. Zarya = 300 Km; Project Magnet = 8 Km for D and I but < .5 Km for F) 
or temporally by determining the mean field at one location over the 
period of several days (e.g. repeat stations) or up to a year (magnetic 
observatories) . Table 3 gives the estimated accuracies used in assigning 
weights. Since the weights for D and I were 'converted 1 to force units 
they should be multiplied by factors of the order of 500 (H or F ~ 30000y) 
for comparison with the other accuracies. Further studies will need to 
be made to ascertain whether these assignments provide a meaningful im- 
provement over unweighted fits. 

4- Results 

Prior to the determination of the final set of coefficients to best 
fit the selected data, preliminary calculations were made to obtain a 
qualitative picture of the meaning of the results and to decide on the 



Table 3 


Estimated accuracies ct for various classes of data 
used in computing the spherical harmonic coefficients 
A. Weights are l/o for H,Z,F, 180 /ctttH for D and 
180 /gttF for I. Dashes appear if elements were not 
observed for this Data Class. 

Elements 


Data Class 

D I 

degrees 

H 

z 

gamma 

F 

Observatory annual means 

.0033 

.006 

5 

15 

15 

Repeat Stations 

.033 

.083 

5 

— 

2 

Land Survey 

.1 

.1 

30 

50 

— 

Shipboard 

.083 

.083 

25 

— 

— 

Ship- towed Magnetometer : 






Proton 

— 

— 

— 

— 

10 

Fluxgate 

_ 

— 

— 

— 

40 

Ai rborne* 

.3 

. 1 

60 

60 

30 

Satellite (Vanguard 3) 

— 

— 

— 

— 

10 


* Observed Elements are: Project Magnet DIF 

Canadian Airborne DHZ 


Univ. of Wisconsin F 
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maximum number of coefficients that would be statistically valid. 

Expanding the coefficients in time about a mean time (~1957) for the 
selected data we may write for gJJ 

<£(t)= g "(t) + (t-i) g“(t) + u-t) 2 g“(t) 

and a similar expression for h™(t). The parameters varied in these 

preliminary calculations were the maximum degree n^. for the spatial 

terms, and n^ for the first and second time derivatives respectively, 

and a parameter k denoting the selection of each kth observation. The 

results that were derived from these fits were the coefficients, their 

first and second time derivatives and the corresponding standard errors 

/ — 2 "' 

o^, a™ and a™ of the amplitudes = V + h^, as previously 

defined. One interesting result from using various non-multiple values 
of k was that for different subsets of the selected data the variations 
in individual parameters did in fact match those estimated by their 
standard errors. This fact gave some confidence that the standard error 
estimates are realistic. The only flaw in this reasoning lies in the 
large gaps in the time - space distributions of the data which allows 
only small variations in the distributions by the choice of various 
values of k. The real test of the validity of the standard error estimates 
can only be made when at least the spatial coverage of data is made 
more uniform. 

The immediate practical problem that arose in the computation of the 
resul Ls was the limitation in the maximum number of coefficients that 
could be used. Including the various check columns, the number of 
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computer core locations required to store the matrix of constants for 
the normal equations was given by N(N-i5) where N = n m (n^ + 2) + 

(n m + 2) + (n^+2) . Thus the storage varied as the fourth power of 

the maximum degree of the coefficients. One alleviating factor that 
was soon apparent was that none of the amplitudes of the second time 
derivatives exceeded the 1.96 a (95$ confidence) level and very few 
attained 0.67 *6 (50$ confidence). This fact resulted in eliminating 
their computation by setting n^=0 for the remainder of the tests. The 
implication here is only that the curvature in the secular change of 
the spatial terms is not significant for the selected data. It is well 
known that at a given location on the earth there is a measurable de- 
viation from linearity over such spans as the 23 years covered by the 
data. It is likely that it would be necessary to take a much larger 
and more uniformly distributed set of data over a longer period before 
a meaningful determination could be made of the second derivatives. 

The final fit to the selected data was thus made with k=l (all data), 
n m =0, n m =5 and n m =7. Since some of the data were anomalous or in error, 
no data were accepted which differed from computed values by more than 
2000/ , a value of the order of five times the standard deviation. Of 
the total number of 21779 component observations this rejection criterion 
only eliminated 84 observations. The resul ting best set of coefficients 
fit to the remaining 21695 observations’, is given in TabLe 4. The 
spatial terms are rounded to the nearest gjimma and the first time 
derivatives to the nearest tenth ‘y/year. Although the average time of 
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Table 4 


n 


m 


8 


g 


h 


1 

1 

2 

2 

2 

3 

3 

3 

3 

4 
4 
4 
4 

4 

5 


7 

7 

7 

7 

7 

7 

7 

7 


0 

-30426 

18.9 


1 

- 2174 

7.3 

5761 

0 

- 1548 

-24.8 


1 

3000 

-0.8 

-1949 

2 

1574 

0.8 

201 

0 

1323 

-0.4 


1 

-2009 

-10.5 

-442 

2 

1275 

3.4 

233 

3 

877 

-1.9 

-118 

0 

957 

0.8 


1 

797 

5.4 

149 

2 

527 

-1.9 

-266 

3 

-400 

-0.2 

-4 

4 

273 

0.8 

-262 

0 

-241 

3.5 


1 

353 

-0.7 

0 

2 

231 

2.5 

124 

3 

-33 

0.6 

-104 

4 

-147 

0.0 

-98 

5 

-79 

1.6 

75 

0 

58 



1 

71 


6 

2 

20 


86 

3 

-241 


58 

4 

-19 


-18 

5 

-0 


-25 

6 

-100 


5 

0 

90 



1 

-47 


-51 

2 

-2 


-22 

3 

-26 


6 

4 

-11 


-38 

5 

26 


44 

6 

6 


-4 

7 

6 


-26 


Best fit to selected data in gamma for the spatial terms 
and y/year for first derivatives. Epoch is 1960.0. 


h 


-1.9 

-14.0 

-17.7 

1.9 
4.0 

-9.0 

-0.9 

-1.7 

3.2 

-5.5 

1.8 

2.9 
- 0.8 
-0.4 
- 0.2 
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the data was 1957.2, the epoch was arbitrarily chosen as 1960.0. The 
values of the amplitudes of the coefficients R and R, and the corre- 
sponding standard errors a and d are given in Table 5 along with their 
ratios R/a and R /a. As can be seen from this table only one of the R 
values and four of the R's fall below the 95$ confidence (1.96c) level. 
From the previous experience in testing the fitting process with lower 
values of n and n^ it is believed that significant coefficients could 
have been obtained for at least one higher degree for both the spatial 
terms and their derivatives if allowed by the computer core size. 

The weighted rms deviations between the field and the data by 
component «re given in Table 6. As in most of these fits it is found 
that the residual to the Z component data is about 50$ higher than for 
the other components. The overall distribution of these deviations 
differs slightly from gaussian in that there are many data beyond 2o 
with over 1$ beyond four standard deviations. This non-gaussian com- 
ponent is of course mainly due to the anomaly 'noise' in the data. 
However, a gaussian curve with a a of about 200y fits the error dis- 
tribution fairly well up to about 400y. 

In order to further investigate how the coefficients fit the 
selected data, weighted averages were taken over 10° blocks in latitude 
and longitude and the resulting signed means printed. The resulting 
numerical 'maps' were qualitatively inspected to determine the areas 
where the computed field was predominantly above or below the measured. 
It was possible to pick out areas of the order of 50° whose deviations 



TABLE 5 


N 

M 

R 

g 


R 

CT 

r /£ 

2 

1 

30426 

9 

3437 

18.9 

1.4 

13.3 

2 

2 

6158 

9 

656 

7.6 

1.3 

5.7 

3 

1 

1548 

8 

203 

24.8 

1.2 

21.4 

3 

2 

3577 

8 

476 

14.1 

1.1 

12.4 

3 

3 

1587 

8 

203 

17.7 

1.1 

16.9 

4 

1 

1323 

6 

204 

•4 

1.0 

•4 

4 

2 

2057 

6 

322 

10.7 

.9 

11.5 

4 

3 

1296 

6 

206 

5.2 

.9 

6.0 

4 

4 

885 

7 

120 

9.2 

.9 

10.7 

5 

1 

957 

6 

169 

0.8 

.9 

.9 

5 

2 

811 

5 

149 

5.4 

.8 

6.6 

5 

3 

591 

6 

103 

2.5 

.8 

3.1 

5 

4 

400 

6 

63 

3.2 

.7 

4.8 

5 

5 

378 

6 

63 

5.5 

.8 

6.8 

6 

1 

241 

5 

49 

3.5 

.6 

5.5 

6 

2 

353 

5 

69 

1.9 

.7 

2.7 

6 

3 

262 

5 

53 

3.8 

.6 

5.9 

6 

4 

109 

5 

22 

1.0 

.6 

1.6 

6 

5 

177 

5 

33 

.4 

.6 

.6 

6 

6 

109 

6 

19 

1.6 

.8 

2.1 

7 

1 

58 

4 

14 




7 

2 

71 

5 

15 




7 

3 

88 

4 

20 





TABLE 5 (Cont.) 


N M 

7 4 

7 5 

7 6 

7 7 

8 1 

8 2 

8 3 

8 4 

8 5 

8 6 


R o R/q 
247 5 49 

26 s 5 

25 5 5 

101 5 20 

99 3 26 

70 4 19 

23 4 6 

27 4 7 

39 4 10 

51 4 13 


8 7 7 4 2 


8 8 




Amplitudes of the.main geomagnetic field R and of their first 
time derivatives R in y and y/year respectively. Standard errors 
o and 6 and the ratio R/o and R/d are also given. The ratios 
should exceed 1.96 for a 95$ confidence and 0.67 for a 50$ con- 
fidence. 


I 


TABLE 6 


RMS deviations by component between computed 
field and selected data with weights based on 
Table 3. 


ELEMENT 

SIGMA (V) 

OBSERVATIONS 

D 

202 

5998 

I 

234 

5087 

H 

211 

2958 

F 

223 

5436 

Z 

344 

2216 

ALL 

214 

21695 
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were systematically a few tens up to a few hundreds of gammas. Thus 
from evidence such as this and the other factors mentioned it is con- 
cluded that there are significant spherical harmonics in the selected 
data of degree beyond seven. 

Illustration of results : 

Since the accuracies of most computed fields are now within a 
percent of the measured values, world maps synthesized from the harmonic 
coefficients all appear very similar. Nevertheless, it is useful to 
attempt to illustrate the scale of the irregularities as represented 
by (.he harmonic descriptions. Maps of F, H, I, and Z for epoch 1965.0 
are given in Fig. 3 as drawn by the techniques previously described by 
Cain and Neilon (1963) • On comparing this Figure with the 1960.0 maps 
given in the earlier paper one notes only slight differences in spite 
of the epoch being five years later and the fact that the present 
charts use 63 instead of 48 spatial coefficients. The reasons for unis 
are that the percentage changes in the field during a five year interval 
are small and that the amplitudes of the additional fifteen coefficients 
are small compared to those of the lower terms. 

Using such a smooth representation it is of course not possible - 
illustrate the numerous small scale anomalies that are found near ihe 
earth's surface. The contours are thus much smoother than are f 
on the standard world magnetic charts which are normally cons'- . a. 
by other methods. One peculiarity of standard charts, however, is 
they am tempt ;,o show as much detail as the data allows with the r~ : 
that there are frequently kinks around anomalies in isolated areas . . 
good data coverage, but smooth contours elsewhere. It is like 




Fig. 3a 





Fig- 3b 






Fig. 3c 



Fig. 3d 
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such maps as shown in Fig. 3 give a more realistic estimate of the 
field averaged over the whole of the earth's surface. 

A prominent feature of the secular change coefficients given in 
Table 4 is the fact that the term does not dominate the series as 
does the g£. Even though there are fewer secular change terms determin- 
able from the data, the resulting secular change maps are at least as 
complex as those for the field itself. The secular change in two of the 
components, F and Z, is drawn in Figure 4. The characteristic feature 
of the F map is that the earth is almost equally divided between positive 
and negative regions. The gross pattern of positive Z is equally simple 
but with the centers of negative change occurring in low latitudes and 
the positives centered only in high latitudes and over Asia. Since the 
estimated confidences in the major secular change terms are far smaller 
than those for the coefficients of the main field (cf . Table 5) the 
details of the patterns illustrated are equally less likely to be signi- 
ficant. 

5 • Comparison With Other Models 

In order to evaluate the way in which the presently derived set of 
coefficients compare with previous models we have made a study of the 
deviations between the available data and the computed field components 
at the points of observation. Complete validation can never be made 
since it is not possible to make comparisons in regions that are void 
of data. Indeed, as new data are accrued it should always be possible 
(provided enough parameters are used) to incorporate this information 





Fig. 4b 
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into the computer programs and improve the field coefficients accordingly! 
The point of this study then can only be to make a comparison of the 
older models to see how they fit the data at various epochs and to see 
whether substantial improvements can be made in the newer models using 
he most recently available data. Less extensive but more detailed 
comparisons have also been done by Heuring (1964) using the set of Van- 
guard 3 data. 

A description of the various models used is given in Table 7. Al- 
though this list does not include all published models, it contains most 
of those widely used in the past with the exception of the 569 coeffi- 
cients derived by Jensen and Whitaker (I960) from the 1955 U. S. Z 
charts. This last model was omitted from the comparison since it has 
been shown in the past to be no substantial improvement over most other 
models in spite of its much larger number of coefficients. 

A statistical comparison was first made between all available survey 
data for the period 1940~1963 and the various models listed in Table 7. 

The data set used was a slightly later version of the one described in 
Tables 1 and 2 so that some data were available for 1963 that were not 
included in the original selection. The rms deviations between the 
observed and computed components were derived as described in the fitting 
of the dala except that no weighting factors were used initially. The 
rms deviations of the angular components AD and AI were computed separately 
in angular units but then weighted by H and F respectively when combined 
with the sums of squares from the AF, AH, and AZ. The data were grouped 


TABLE 7 


Symbol 

A 

LME 


F+L 

F+L(t) 


J+C 

N+0 


V 

J+M 

USSR 


F+K 


Reference* 

Set of coefficients listed in Table 4 
Leaton, Mai in and Evans . 1964* (A fit 
with = 8 and ^ = 7 to a set of X and 
Y computed data assuming a spherical earth) . 
( Finch and Leaton . 1957) 

The (t) added to F+L indicates that the 
secular change terms of Leaton (1962) were 
used. 

Jensen and Cain (1962) 

( Nagata and Oguti . 1962) The use of the (2) 
means that both the internal and external 
fields were added. If (i) only appears, the 
evaluation was made using only the internal 
field. 

Vestine (i960) 

Jones and Melotte (1953) 

Adam et al (1962) Table 4, Adam et al (1963) 
secular change from appendix (61 observatories 
1954-1959) 

Fauselau and Kautzleben (1956) 


*A listing of some of these coefficients can be found in a report 
by the authors ( Cain et al . 1964) 
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according to each year 1940-1963 and by the source of the observations. 

The source categories were: observatory annual means, surface (land or 
sea surveys), airborne, and satellite (Vanguard-3 F only). Results 
were tabulated for each of the component-source-year groups. Summaries 
were also tabulated by source-year and year alone and also for each 
source-component or source alone for all years 1940-1963. Using the 
set of coefficients A, data were first rejected if their deviation 
exceeded 3000y in the components, F, H, Z or 10° in the angles I or D. 

In each succeeding comparison the same rejection criteria were maintained 
for the remaining data. A summary of the results by year is given in 
Table 8 for the eleven different models. The last column is labeled 
"Maximum number of observations" since in order to conserve computer 
Lime, only a fraction of the available data were used for each model. 

The figures in the table represent the rms deviations from the largest 
data sample in those cases where different sample sizes were used. In 
making such comparisons using different subsets of the data the results 
varied up to a few tens of gamma. Thus only those differences exceeding 
this amount should be considered significant. 

The two set of field coefficients, A and LME, appear to have the 
lowest residuals to the data of all coefficients sets over most of the 
interval 1940-1963. This improvement is striking even over the F-t-L(t) 
coefficients which are the next best set over the whole interval 1940-1964. 
Particular other sets of coefficients show relatively low residuals 
near their epoch time but not elsewhere. Noteworthy in this respect are 
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Vestine 1 s (I960) 1945-epoch coefficients which appear to give consistently 
low residuals over the years 1942-1945* The N+0(i) set appear to give 
consistently good results in the interval 1952-1956 and for a few pre- 
vious years but tend to become worse at and after their epoch of 1958.5. 
This tendency for field models to better fit the data prior to their 
epoch is the result of the fact that most authors are attempting to 
compute a field valid at the time of publication based only on past data. 
This result indicates that 'correcting' data to a given epoch is not 
on the whole successful. 

The set of coefficients derived by Nagata and Qguti (1962) were used 
to determine how the addition of external terms change the residuals to 
the data. As can be seen from the column N+0(2) and N+0(i) , the use of 
internal terms only gives a better fit. The reason for this discrepancy 
is probably that the external terms represent only the errors in the 
maps from which the determinations were made and are unrelated to any 
realistic external field. 

The comparison between the Finch and Leaton (1957) coefficients with 
and without Leaton' s (1963) time derivatives is indicative that the use 
of secular change terms is important even over a few years. The in tent 
of the computation made by Jensen and Cain (1962) was to improve the 
residuals over those of F+L (no time derivatives) . It can be seen that 
for the years 1959-1963 the J+C residuals are indeed less that the F+L 
by about 10 %. However, when the time derivatives are added (F+L(t)), this 
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improvement is eliminated so that the two sets of coefficients are 
merely comparable . 

Having established th..t the presently derived coefficients A and 
LME clearly give results that are much more representative of the field 
than any of the others it is useful to compare these in more detail to 
see where the differences may lie. In the yearly comparisons of Table 
8 the residuals are smaller for the A field than LME for only nine of 
the 24 years and except for the 1 9 6 i data are all grouped near the 
beginning of the interval. All of these comparisons were made using the 
survey data with equal weights whereas the fit resulting in the model A 
was actually made with weights established by Table 3. Table 8 was con- 
structed in this unweighted way in order to simplify the discussion and to 
see whether the use of weights made any appreciable difference in the 
apparent relative merits of the various fields. If the data are instead 
weighted according to Table 3 and a comparison again made of the residuals 
against the various computed fields, the A and LME sets again have the 
lowest residuals by a factor of the order of two from the others. The 
relation between A and LME however is altered in that the yearly residuals 
for A nearly all become the smaller. As summarized in Table 9 the difference 
averages a third for 1940-1955 but only 7% for the data after 1955. It is 
thus apparent that at this point the relative merits of the two sets of 
coefficients depend on whether the weights are accepted as realistic. 

In both Tables 8 and 9 the data for 1940 are fit much better by A 



TABLE 9 


A comparison of the rms residuals between observed data and data 
synthesized using fields A and LME and the weighting factors based 
on Table 3. 



rms 

Residuals 

(gamma) 



Years 

A 

LME 

Years 

A 

LME 

1940 

286 

433 

1952 

217 

244 

1941 

237 

336 

1953 

204 

207 

1942 

225 

278 

1954 

245 

242 

1943 

241 

323 

1955 

197 

206 

1944 

240 

260 

1956 

232 

223 

1945 

209 

281 

1957 

222 

256 

1946 

243 

277 

1958 

258 

272 

1947 

201 

254 

1959 

140 

154 

1948 

218 

253 

I960 

249 

270 

1949 

183 

217 

1961 

225 

236 

1950 

290 

307 

1962 

264 

258 

1951 

246 

242 

1963 

239 

266 

1940-1955 

249 

322 




1955-1963 

210 

224 
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than by LME. The majority of data for this year is that from the Soviet 
Union. As there was only observatory data since this time from that 
area there is possibly a major difference between A and LME over Asia. 
Since Lea ton . Mai in and Evans (1964) indicate that they used for their 
analysis the data taken from the USSR 1950 T(=F) and 1955 D, H and Z 
charts, one may conclude that there are some' large differences between 
those charts and the total set of digital data on file. One other 
possibility is that the use of only single time derivatives may not be 
adequate for representing the data for this region weighted as heavily 
as it is with data at 1940. 

Table 9 indicates fairly low residuals for weighted data for both 
models for 1959- The reason for this is the inclusion of the Vanguard 
3 data which have relatively high weights and match either model to ~ 50 y 
rms. The fact that it is possible to fit these satellite data to a 
precision of 17y ( Cain and Hendricks , 1964) is further indication that 
higher order harmonics are needed. One possibly significant difference 
between the A and LME models is that whereas the rms deviations to the 
Vanguard 3 data are almost the same, the mean errors for A are close to 
zero whereas the measured F average about 20y below that predicted by 
LME at all altitudes. We cannot offer any clear physical or mathematical 
reason for this discrepancy except to reiterate some of the differences 
in derivation that might be significant. The most apparent difference 
of course is that the fit to derive A included the Vanguard 3 data 
whereas LME did not. Thus recalculating the field in the volume of the 
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Vanguard 3 measurements becomes an extrapolation for LME but an inter- 
polation using tfye A coefficients. Secondly, the LME field was based 
on the spherical earth approximation which we have already shown can 
give rise to errors at altitude. 

Bearing in mind this average difference between A and LME using 
one sample of low altitude satellite data one may also consider how 
the field models extrapolate to larger distances. One such sample of 
data on which such comparisons can be made is the small amount of 
Explorer X total field data available for March 25, 1961 from 1544 to 
1740 U.T. which range from dipole latitudes (cf. Chapman . 1963) -8 to 
-31° and altitudes from 1.7 to 7 earth radii. As can be seen in Fig. 

22 of the paper by Heppner et al (1963) the measured field was below 
that computed by F+L by some 40y at the lowest altitude but became equal 
to and greater than that computed as the altitude was increased. Since 
the lowest altitude portion of the curve will show the least distortion 
from the effects of trapped particles and the magnetopause we use it 
to compare with the field models. The maximum AF (measured minus com- 
puted) for seven of the field models is given in Table 10. It is in- 
teresting to note that the AF for all four fields LME, J+C, N+0(i), and 
F+L(t) is very nearly -30y whereas the lowest deviation of -lBy occurs 
for A. Ignoring the altitude uncertainty mentioned by Heppner et al 
(1963) for the measurement it would appear that the A model gives the 
best prediction of the main field. 

Another independent set of information relating to the way in which 
various field models fit a low latitude field line is provided by 



Table 10 > 


Difference AF (measured minus computed) between total field 
observed by Explorer X at 14°S0UTH, 4°WEST, 5000 Km altitude, 
arv:: field c c“\:-ed by several models . 


F ield 

AF 

A 

-18 

LME 

-30 

J+C 

-30 

N+0( i) 

-29 

F+L 

-42 

F+L( t) 

-30 

USSR 

-92 


Table 11 


Distance R (km) between position of observed artificial aurora 
( Leonard. 1963) and traces to 100 Km altitude using field models indicated. 


Model 

R(Km) 

F+L 

56 

J+C 

46 

LME 

35 

A 

26 
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Leonard 1 s (1963) observations of the artificial aurora conjugate to one 
Johnson Island bomb tests. In his Figure 3 Leonard shows the position 
of the observed aurora and the E-layer intersections traced with various 
field models. Assuming an approximate injection point of 16.50°N, 
169.64°W and 400 Km altitude (which appears to match the positions given 
in his Fig. 3) we tabulate in Table 11 the distance from the observed 
position given by Leonard (l6.6°S, 175»82°W) to the point intersecting 
100 Km altitude as traced from the estimated injection point, using 
four field models. As can be seen in this table, the A trace matches 
the observed position with the smallest error. 

In order to obtain a detailed picture of the way in which the models 
fit the near surface data we have plotted some of the deviations for a 
few of the ship and aircraft data. A map showing the tracks illustrated 
is given in Figure 5. Most of these were chosen to be near the regions 
void of data and should thus result in the largest deviations. The de- 
tailed plots of the differences of the observations from the A and LME 
models for each of these tracks are given in Figure 6. The differences 
shown are the measured less computed as observed except for the angular 
observations AD and AI which were again converted to force units. There 
are several useful conclusions that can be gleaned from such plots in 
regard to the way in which these two model fields fit the observations. 
The most striking feature is that there are long-period systematic 
deviations for either field model which indicate that either more or 
better coefficients need to be derived. The rms deviations given in 
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Table 9 thus appear to reflect significant deviations over large distances 
in addition to the scatter in the observations caused by short-wavelength 
surface anomalies. Although there are many instances where the data 
deviations from the two different models are merely different, there 
are more instances of similarities. For example, for MAGNET ( 316) the 
computed D is too low by about a degree for both A and LME. One also 
finds that the structure of the deviations is often quite similar and 
that the wavelengths are less than 360/n^ so as to again confirm that 
higher order harmonics are indeed necessary. Such examples can clearly 
be seen in MAGNET (633) AF and FAI and most of the ZARYA plots. 

One can also compare the secular change observed at the magnetic 
observatories with that predicted by the A and LME models. Sample plots 
of the annual mean observations at three observatories Alibag, San Juan, 

and Sodankyla are given in Fig. 7. The dots represent the observations, 

the solid lines the field components predicted by A and the dashed lines 
those predicted by LME. The LME model was based on data estimating the 
secular change at 1965.0 whereas the A model used data selected over the 
whole interval 1940-1962 but weighted towards more recent times so that 
the mean date was 1957.2. Thus one would expect the slopes of LME 
computed values to match those at 1965*0 while the A values should 
approximate the mean slope over at least the decade of the 1950's. As 

seen in Fig. 7 this tendency is evident only for the H and D curves at 

Sodankyla and Alibag. The absolute differences between the annual means 
at the three observatories and that computed by the two field models 
are of the same order as one might expect at three surface positions. 
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Conclusions 

We have presented a new model of the earth 1 s main magnetic field 
and secular change which is a significant improvement over all past 
models for the period 1940-1963 with the exception of one recently derived 
by Lea Lon t Mai in , and Evans (1964) • This last model shows a comparable 
fit lo Lhe data but does not extrapolate as well with altitude perhaps 
partly due to having used Lhe spherical approximation to the earth* s 
surface. The allowance for the oblateness of the earth was illustrated 
to be a necessary refinement at Lhe levels of accuracy now being reached. 

It was apparent in comparing both of the two models against survey 
data that there are significant higher order harmonics in the field. 

The deviations from the fits appear to arise at least as much from long 
wavelength systematic variations as from the small-scale anomaly n noise n . 
Prior conclusions that coefficients of degree exceeding n-6 are not 
meaningful ( Fanselau and Kautzleben , 1956) appear to result from analyzing 
charts instead of the data themselves. We must conclude that a better 
field model can be derived by analysis of the raw data than any existing 
world chart. It appears that better maps could be derived from such 
analyses as shown here than from more conventional techniques whereby 
the data are "corrected" to a given epoch and field contours drawn 
manually for each component. 

The present model does not allow for external fields nor did the 
evaluation make any corrections to the data for magnetic, disturbance. 

A future calculation is planned to see whether any meaningful improvements 
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in the model could be made by the use 
for the systematic depressions in the 
magnetic disturbance ( Sugiura , 1964) • 
also include the use of higher degree 


of external terms and adjustments 
field during and following 
These further improvements will 
terms as necessary to match the 


complexities of the field. 
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